s 


NPS-MA-93-018 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


SATELLITE  MOTION  AROUND  AN  OBLATE 
PLANET:   A  PERTURBATION  SOLUTION  FOR 

ALL  ORBITAL  PARAMETERS 

by 

D.  A.  Daniel  son 

G.E.  Latta 

C.P.  Sagovac 

S.D.  Krambeck 

J.R.  Snider 

Technical  Report  For  Period 
January  1993  -  April  1993 


Approved  for  public  release;  distribution  unlimited 
Prepared  for:   Naval  Postgraduate  School 
Monterey,  CA  93943 


FedDocs 

D  208.14/2 

NPS-MA-93-018 


^"Voes-^-^-^ 


NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  93943 

Rear  Admiral  T.A.  Mercer  Harrison  Shull 

Superintendent  Provost 

This  report  was  prepared  in  conjunction  with  research  conducted  for  the  Naval  Postgraduate 
School  and  funded  by  the  Naval  Postgraduate  School. 

Reproduction  of  all  or  part  of  this  report  is  authorized. 

This  report  was  prepared  by: 


Ha<:<;ifip=d 


TY  CLASSIFICATION  OF   THIS   PAGE 


REPORT  DOCUMENTATION  PAGE 


form  Approved 
OMB  No  0704  0188 


— 


9ORT  SECURITY  CLASSIFICATION 

classified 


lb    RESTRICTIVE  MARKINGS 


":■•■  O 

:;"  ■:'■, 

:  | 

O  r 

w  Z 

d  o 

C  ?  D  i 

■,  5 

1  n 

m  < 

w 

o 

X 

O 

o 

i- 

1URITY  CLASSIFICATION  AUTHORITY 


CLASSIFICATION  /  DOWNGRADING  SCHEDULE 


3     DISTRIBUTION /AVAILABILITY  OF   REPORT 

Approved  for  public  release; 
distribution  unlimited 


— 


_ 


ORMING  ORGANIZATION  REPORT  NUMBER(S) 


S-MA-93-018 


5    MONITORING  ORGANIZATION  REPORT  NUMBER(S) 

NPS-MA-93-018 


ME  OF  PERFORMING  ORGANIZATION 


val  Postgraduate  School 


6b    OFFICE   SYMBOL 
(If  applicable) 

MA 


la    NAME  OF   MONITORING  ORGANIZATION 


Naval  Postgraduate  School 


)RESS  {City,  State,  and  ZIP  Code) 

mterey,  CA  93943 


7b    ADDRESS  (City,  Sfafe    and  7IP  Code) 

Monterey,  CA  93943 


ME  OF  FUNDING /SPONSORING 
5ANIZATION 

/al  Postgraduate  School 


8b    OFFICE    SYMBOL 
(If  applicable) 

MA 


9    PROCUREMENT  INSTRUMENT   IDENTIFICATION  NUMBER 


OM,N 


5RESS  (City,  State,  and  ZIP  Code) 

nterey,  CA  93943 


10    SOURCE    Of    FUNDING  NUMRt  RS 


PROGRAM 
ELEMENT   NO 


PROJECT 
NO 


TASK 
NO 


WOPf    UNIT 
ACCESSION  NO 


E  (Include  Security  Classification) 

sllite  Motion  Around  an  Oblate  Planet:  A  Perturbation  Solution  for  all  Orbital  Parameters 


£orr3ameK8n,S)G.E.  Latta,  C.P.  Sagovac,  S.D.  Krambeck,  J.R.  Snider 


PE  OF  REPORT 

hnical 


13b    TIME   COVERED 

from  1-93       to  4-93 


14    DATE   OF   REPORT    (Year,  Month,  Day) 


6  May  93 


IS    PAGE   COUNT 


24 


PLEMENTARY  NOTATION 


COSAii  CODES 


LD 


GROUP 


SUB  GROUP 


18    SUBJECT    TERMS  (Continue  on  reverse  if  necessary  and  identity  by  block  number) 

Perturbation  solution,  oblate  planet,  orbital  parameters 


TRACT  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

The  search  for  a  universal  solution  of  the  equations  of  motion  for  a  satellite  orbiting  an  oblate 
let  is  a  subject  that  has  merited  great  interest  because  of  its  theoretical  and  practical  implications, 
e,  a  complete  first-order  perturbation  solution,  including  the  effects  of  the  J2  terms  in  the  planet's 
sntial,  is  given  in  terms  of  standard  orbital  parameters.   The  simple  formulas  provide  a  fast  method 
predicting  satellite  orbits  that  is  more  accurate  than  the  two-body  formulas.   These  predictions  are 
wn  to  agree  well  with  those  of  a  completely  numerical  code  and  with  actual  satellite  data.   Also,  in 
appendix,  it  is  rigorously  proven  that  a  satellite  having  negative  mechanical  energy  remains  for  all 
i  within  a  spherical  annulus  with  radii  approximately  equal  to  the  perigee  and  apogee  of  its  initial 
llating  ellipse. 


RIBUTION  /AVAILABILITY  OF   ABSTRACT 
MCLASSIFIED/UNLIMITED       D  SAME   AS  RPT  □  DTIC  USERS 


21     A 


BT3tfc7assfff&f 


CLASSIFICATION 


■ME  OF  RESPONSlBlE   INDIVIDUAL 

i.  Daniel  son 


?2b  TELEPHONE  (Include  AieaCode) 

408-656-2622 


?2c    OFFICE   SYMBOL 

MA/Dd 


*n1473,  JUN  86 


Previous  editions  are  obsolete 
S/N   0102-LF-014-6603 


SECURITY    CLASSIFICATION  OF    THIS   PAGE 


SATELLITE  MOTION  AROUND  AN  OBLATE 

PLANET:  A  PERTURBATION  SOLUTION  FOR 

ALL  ORBITAL  PARAMETERS 

D.  A.  Danielson,  Professor 

G.  E.  Latta.  Professor 

C.  P.  Sagovac,  LT,  U.S.  Navy 

S.  D.  Krambeck,  LT,  U.S.  Navy 

J.  R.  Snider,  LTC,  U.S.  Army 

Mathematics  Department 

Naval  Postgraduate  School 

Monterey,  California  93943 

Abstract 

The  search  for  a  universal  solution  of  the  equations  of  motion  for  a 
satellite  orbiting  an  oblate  planet  is  a  subject  that  has  merited  great 
interest  because  of  its  theoretical  and  practical  implications.  Here,  a 
complete  first-order  perturbation  solution,  including  the  effects  of  the 
J2  terms  in  the  planet's  potential,  is  given  in  terms  of  standard  orbital 
parameters.  The  simple  formulas  provide  a  fast  method  for  predict- 
ing satellite  orbits  that  is  more  accurate  than  the  two-body  formulas. 
These  predictions  are  shown  to  agree  well  with  those  of  a  completely 
numerical  code  and  with  actual  satellite  data.  Also,  in  an  appendix,  it 
is  rigorously  proven  that  a  satellite  having  negative  mechanical  energy 
remains  for  all  time  within  a  spherical  annulus  with  radii  approximately 
equal  to  the  perigee  and  apogee  of  its  initial  osculating  ellipse. 


1     Introduction 

A  characteristic  feature  of  practical  orbit  prediction  is  that  the  engineer  may  deal  with 
numerous  satellites  in  a  great  variety  of  orbits.  Under  these  circumstances  analytical  relations 
which  can  quickly  approximate  an  orbit  may  be  far  superior  to  large  numerical  programs. 
While  many  analytical  models  have  been  developed  for  the  artificial  satellite  age,  most  are 
not  used  in  practical  orbit  prediction  because  they  violate  one  or  more  of  the  following 
principles: 

•  The  method  should  provide  a  solution  that  is  significantly  more  accurate  than  the 
two-body  solution. 

•  The  real  physical  effects  of  the  orbit  should  be  easily  distinguishable  in  the  solution. 

•  The  solution  should  be  universal:  it  should  be  valid  for  all  orbital  parameters. 

The  problem  of  predicting  the  motion  of  a  satellite  perturbed  only  by  the  oblateness  of  the 
planet  has  received  considerable  attention  following  the  first  launchings  of  artificial  satellites 
about  the  Earth.  Some  of  the  studies  of  this  problem  by  means  of  general  perturbation 
theories  are  listed  at  the  end  of  this  paper.  Techniques  have  involved  expansions  in  powers 
of  y/Jl-  averaging  processes,  the  use  of  spheroidal  coordinates,  and  the  edifice  of  Hamiltonian 
mechanics.  It  is  not  the  intention  of  this  present  paper  to  compare  the  various  methodologies 
used.  Suffice  it  to  say  that  many  researchers  believe  a  solution  which  embodies  all  of  the 
above  principles  was  not  achieved  (e.g.,  see  Taff). 

The  basic  procedure  used  in  this  paper  to  solve  the  differential  equations  of  motion  is 
the  perturbation  technique  known  as  the  Method  of  Strained  Coordinates.  This  technique 
was  first  applied  to  the  title  problem  by  Brenner.  Latta,  and  Weisfield.  Using  a  mean  orbital 
plane  to  specify  an  arbitrary  orbit,  they  were  only  able  to  obtain  a  partial  solution  (e.g.,  the 
eccentricity  was  assumed  small  and  initial  conditions  were  not  considered). 


Here  we  use  coordinates  in  the  true  orbital  plane  to  cast  the  differential  equations  into  a 
simplified  form,  as  was  originally  done  by  Struble. 

2     Orbital  Kinematics 

Figure  1  shows  the  usual  reference  system  of  spherical  coordinates  (r,  a,  0).  The  radial 
distance  r  is  measured  from  the  center  of  the  planet  0  to  the  satellite  S.  The  line  O7  is  in 
a  direction  fixed  with  respect  to  an  inertial  coordinate  system.  The  right  ascension  a  is  the 
angle  measured  in  the  planet's  equatorial  plane  eastward  from  the  line  O7.  The  declination 
or  latitude  (5  is  the  angle  measured  northward  from  the  equator.  The  position  vector  r  of 
the  satellite  in  the  spherical  coordinate  system  is 

r  =  r(cosQ  cos  $)hi  +  r(sin  a  cos/?)b2  +  r(sin  l3)b3  (1) 

where  (bi,b2,b3)  are  orthonormal  base  vectors  fixed  in  the  directions  shown. 

We  can  also  locate  the  satellite  by  its  polar  coordinates  (r.  6)  within  a  (possibly  rotating) 
orbital  plane  that  instantaneously  contains  its  position  and  velocity  vectors.  Here  6  is  the 
argument  of  latitude,  i.e.,  the  angle  measured  in  the  orbital  plane  from  the  ascending  node  to 
the  satellite.  The  orbital  plane  is  inclined  at  an  angle  i  to  the  equatorial  plane  and  intersects 
the  equatorial  plane  in  the  line  of  nodes,  making  an  angle  Q  with  the  0~)  line. 

We  introduce  another  orthonormal  set  of  base  vectors  (Bi.B2,B3)  which  move  with  the 
satellite  so  that  Bi  is  in  the  direction  of  the  position  vector  r.  B2  is  also  in  the  orbital  plane. 
and  B3  =  Bi  x  B2.  The  basis  (bi,  b2.  b$)  may  be  transformed  into  the  basis  (Bj.  B2.  B3)  by 
a  succession  of  three  rotations.  First  the  basis  (b!,b2.b3)  is  rotated  about  the  b3  direction 
by  the  angle  ft,  next  the  basis  is  rotated  about  the  new  1-direction  by  the  angle  i,  and 
finally  the  basis  is  again  rotated  about  the  new  3-direction  by  the  angle  6.  The  two  sets  of 
base  vectors  are  related  by  the  product  of  the  rotation  matrices  representing  each  successive 


rotation  (as  explained  in  the  book  by  Danielson^ 


B2 
B3 


cos  6    sin#    0 

—  sin  0    cos  6    0 

0  0       1 


1         0  0 

0        cos  i     sin  i 
0     —  sin  i     cos  i 


cos  n   sin  n   0 

sin  ft     cos  ft     0 
0  0        1 


b2 
b3 


or 


B, 
B2 
B3 


cos  6  cos  ft  —  sin  6  cos  i sin  ft        cos  0  sin  ft -f  sin  #  cos  ?  cos  ft     sin  0  sin  ? 
—  sin  6  cos  ft  —  cos  6  cos  i  sin  ft     —  sin  9  sin  ft  4-  cos  6  cos  z  cos  ft    cos  6  sin  i 
sin  i  sin  ft  —  sin  z  cos  ft  cos  i 


(2) 


bi 
b2 
b3 


The  position  vector  r  has  only  one  component  in  the  rotating  basis: 

r  =  rB,  (3) 

Using  the  first  of  equations  (2),  we  obtain  the  components  of  r  in  the  fixed  basis: 

r  =  r(cos#cos  ft  —  sin  6  cos  i  sin  ft)^ 

+  r(cos#sin  ft  +  sin  6  cos  ?  cos  ft)b2  +  r(sin  ^sin  ?)b3  (4) 

Equating  the  components  of  equations  (1)  and  (4).  we  can  obtain  the  following  relations 
among  the  angles  (a,  j3)  of  the  spherical  coordinate  system  and  the  astronomical  angles 
(i,ft,0): 

sin  3  =  sin  #sin  ?  (5) 

cos  3  —  cos  6  sec( q  —  ft ) 

The  velocity  dr/dt  of  the  satellite  is  obtained  by  differentiating  (3)  with  respect  to  the 

time  t: 

dr         dr  JR, 

(6) 


dr       dr  rfBj 

-di  =  TtBi  +  r^r 


Since  the  orbital  plane  must  contain  the  velocity  vector,  we  have  to  enforce 

dB, 


dt 


B,  =  0 


(7) 


Substitution  of  equations  (2)  into  equation  (7)  leads  to  a  relationship  which  uncouples  the 

equations  for  ft(0)  and  i{6): 

dQ        tantf    di 

(8) 


dO        sin?    dO 


The  velocity  (6)  can  then  be  written 


dv       <fr  d6  „        .di\n 

—  =  — B!  +  r—     1  +  tan<9cot?  —    B2 
dt       dt  dt  \  dO 


(9) 


In  the  following  part  of  this  paper,  we  will  obtain  expressions  for  r(0),  i{9),  fi(0),  and 
dt/dO(0).  The  position  and  velocity  vectors  of  the  satellite  then  may  be  calculated  from 
the  formulas  in  this  section.  The  classical  orbital  elements  p,  e,  and  u;  are  the  semilatus 
rectum,  eccentricity,  and  argument  of  perigee  of  the  instantaneous  (osculating)  conic  section 
determined  by  the  position  and  velocity  vectors.  If  needed,  p(9),  e(0),  and  u>(9)  can  be 
obtained  from  our  solution  r(9)  and  dt/d9(9): 


V- 


GM  {%) 


t  cos(9  —  u:)  = 1 

r 

esin(»-u.-)  =  ^(J 


3     Equations  of  Motion 


The  expressions  for  the  kinetic  and  potential  energies  per  unit  mass  of  a  satellite  orbiting 
around  an  oblate  planet  are  respectively: 

2  /,i\2  /,\2" 


T  =  \ 


%)+*(%)+<>«*>$ 


V  =  - 


GM 


1  + 


J2R2 

2r2 


(l  -3sin2/3) 


(10) 


11) 


where  G  is  the  gravitational  constant,  M  is  the  mass  of  the  planet,  R  is  the  equatorial  radius 

of  the  planet,  and  J2  is  the  constant  coefficient  of  the  spherical  harmonic  of  degree  2  and 

order  0  in  the  planet's  gravitational  field.    Substitution  of  these  equations  into  Lagrange's 

equations 

d  9(7- -V)       d 

S^Ff-^(r-,)  =  °  1  =  r,a,  or  0 


results  in  the  following  equations  of  motion: 

d>r         (d3\2  2,/^V         9Y 

__r(_j    _rcos^^-j    =  -  —  (12) 

i^s^-^dr-s       <i3) 

Initial  conditions  are  established  by  requiring  that  at  the  initial  time  to  the  orbital  pa- 
rameters of  the  usual  two-body  orbit,  the  conic  section  determined  by  the  initial  position  and 
velocity  vectors,  are  known.  The  actual  orbit  is  then  tangent  to  this  initial  instantaneous 
conic  section  at  to  (see  Figure  1).  Equating  the  initial  position  and  velocity  vectors  given  by 
equations  (3)  and  (9)  to  the  two-body  expressions,  we  obtain 

r«o)  =  — ^ :,  (14) 

1  +  e0cos(#0  -  u;0) 

dr  eohos\n{0o  -  uj0)  ,... 

-r;(to)  = (Id 

at  p0 

dO  h0 

-«o)  =  Tr -1— -  (16) 


dt->       rl 


1  +  tan0ocoUog(0o; 


Wo)  =  io  (17) 

n(^0)  =  o0  (is) 

Here  h0  =  yJGMpo  is  the  initial  value  of  the  satellite's  specific  angular  momentum  about  the 
center  of  the  planet,  and  the  subscript  0  on  a  symbol  denotes  that  the  parameter  is  evaluated 
at  the  initial  time  t0. 

We  immediately  have  two  integrals  of  the  equations  of  motion: 

T  +  V  =  constant  (19) 

rz  cos^  8—  =  constant  (20) 

dt  v      ; 


Equation  (19)  simply  states  that  the  mechanical  energy  of  the  satellite  remains  constant. 
Now,  from  equations  (1)  and  (16) 


.2  „_2 


da 


dr 


r  cos   0—  =  r  x  —  •  b3  =  h0  cos  z0 


(21) 


dt  '  dt 

Equation  (21)  simply  states  that  the  component  along  the  polar  axis  of  the  specific  angular 

momentum  of  the  satellite  remains  constant.  Inserting  equations  (3)  and  (9)  into  equation 

(21),  we  obtain 

di  r    m<s  if  di\ 

(22) 


dt         r  cos  if  dx 

—  = 1  +  tan  6  cot  z  — 

dv       h0  cos  z0  \  dB 


This  allows  the  independent  variable  to  be  changed  from  t  to  6. 

Letting  u  =  po/r,  and  using  equations  (5).  (21).  and  (22).  we  can  rewrite  the  remaining 
equations  of  motion  ( 1 2)— ( 13): 

di        —  2J  u  sin  6  cos  8  sin  i cos2  i 


d6  '        -^-  +  2Jus\n26cos3i 


COS  I 


(23; 


d2u  cos2  ?       J  cos* 

+  w  =  — -r-  + r 


d0< 


c* 

d2u 


2  :  r  fa 

u2(\  -  3sin20sin2?)  +  2u  —  sin  0  cos  0(1  -3  cos2?; 


du 


Au  — —  sin2  6  cos2  z '  —  2  I  —  )    sin2  6  cos2  z 


d6 
AJ2us'm3  6cos6i 


du  .  2  d   f    du\  2  ■ 

u  —  cos  6/(2  +  sin   z)  +  —  I  u  —  I  sin  6*  cos   ? 

The  terms  in  (24)  with  d2u/d92  can  be  combined,  yielding  the  equivalent  equation 
d2u  [  cos2  i       J  cos2  i .  0 ,  .  0  „  9  . 

^7  +  u  =  \  —  +  — — \-u  (]  + sin  ^7 cos  -  - 3)) 

du  o  /  c/u  \         o  t  , 

+2u  —  sin  0  cos  0(1  -3  cos2?)  -2    —      sin2  0  cos2?] 
dv  \  du  I 


+ 


4J2izsin30cos6?' 


22  du  2  I  du  \  2 

u   sin p cos   i '  —  u—  cos 6(2  +  sin   z)  —  I  —  J    sin 6 cos   z 
u0  \  d6  I 


f        4 J??  sin2  6  cos4  ?       AJ2u2  sin4  0  cos8  ? 
-5-1  + + : 


[24) 


(25) 


Here  we  have  introduced  the  shorthand  notation  c  =  cos?o,  5  =  sinz0.  J  =  SJ2R2 I~Pq- 

4     Perturbation  Procedure 

The  differential  equations  (23)-(24)  are  coupled  by  the  nonlinear  terms  and  apparently 
cannot  be  solved  analytically.  If  we  expand  the  right  sides  of  (23)  and  (25)  in  a  Taylor  series 
expansion  in  powers  of  J  and  retain  only  terms  up  to  order  J2,  the  equations  simplify  to 

di        —2Jus'm6cos6s\nicos3i      4  J2u2  sin  i  cos7  i    .   ,  -        „       _,.  To.  .__> 

—  =  = + sin3  9  cos  6  +  0(J3)  (26) 

dv  cl  c4 

d2u  cos2 1       J  cos2  7  ( —Aus'm2  6  cos4  i         9r  .   9  „.„       ,.      „., 

m  +  u  =  —  +  — r-{ 5 +  «  P  + sin  '( J  cos  '  - 3)] 

+  2u-^sin0cos6»(l  -  3  cos2  ?)  -2  ( -^  ]    sin2flcos2?)  (27) 

do  \  dv  J  ) 

4J2us\n2  Ocos6 1  (   9r  „   .   ,  „  „       ,  ...      3u  sin2  #cos4  i 

+ \u2[-\  +3sin20(l  -2cos2?)  + 

c4  I  cz 

+  U-377  sin  0  cos  0[7  cos2?  -  5]  +  (  ~  |    sin20cos2?i  +  0{J3) 
dv  \  dv J  J 

Here  the  term  in  the  0  symbols  indicates  that,  for  all  sufficiently  small  J,  the  error  is  less 
than  a  constant  times  J3.  The  equations  (26)— (27)  are  identical  to  those  used  as  the  starting 
point  in  the  analysis  of  Eckstein,  et  al. 

It  is  reasonable  to  expect  that  the  solution  for  u  will  be  arbitrarily  close  to  the  two  body 
solution.  1  +  tQCOs(6  —  u;0),  when  J  is  close  to  zero.  This  assumption  is  consistent  with 
letting 

u  =  \  +  e0cosy  +  Jui  +  J2u2  + ...  (28) 

y  =6-u0  +  Jyi  +  J2y2  +  ...  (29) 

i  =  i0  +  JU  +  J2*2  +  •  •  •  (30) 

An  algorithm  for  the  perturbation  procedure  is: 

Let  n  =  \ 

Substitute  expressions  (2b)-(30)  into  the  equations  of  motion  (26)-(27) 


Equate  the  coefficients  of  Jn 

Choose  the  arbitrary  constants  so  secular  terms  will  not  arise. 

Solve  for  the  nth  order  solution 

Satisfy  the  initial  conditions  (lj)-(18) 

Iterate  on  n 

The  calculations  were  carried  out  with  the  symbolic  manipulation  program  MACSYMA. 
In  this  paper  we  only  briefly  outline  these  calculations;  for  more  details  see  the  theses  of 
Sagovac  and  Snider. 

Beginning  by  substituting  equations  (28)  and  (30)  into  (26),  and  equating  the  terms 
multiplied  by  J,  we  obtain 

J-  =  -scsm20  -  ^p  sm(y  +  20)  +  ^sm(y  -  20)  (31) 

A  solution  to  this  equation  is 

i,  =  ^  cos 20  +  ^cos(y  +  20)  +  ^cos(y  -  20)  +  h\  cos(2y  -  20)  +  K2  (32) 

2d  2 

The  last  two  terms  may  be  added  because  they  are  to  lowest  order  homogenous  solutions 
to  equation  (30).  The  term  multiplied  by  the  constant  K\  was  added  to  eliminate  secular 
terms  in  i2\  note  that  differentiating  this  term  with  respect  to  0  produces  terms  multiplied 
by  J,  from  equation  (29).  The  constant  K2  was  added  to  satisfy  the  initial  condition  (17). 
which  implies  that  i\(0o)  =  0  so 

K2  =  -—  cos  20Q —  cos(30o  -  w0)  -  "—r-  cos(0o  +  u;0)  -  A'i  cos  2cc?o 

2o  2 

Substituting  equations  (28)— (30)  and  (32)  into  (27).  and  equating  terms  multiplied  by  J 

yields 

tPu,  3.s2        „  /    ,^2        \       1 
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cPu-i  3s2        ,  /    5s2        \       1 .  oo         o, 

-Jgt  +  «i  =  1  "  ~Y  +  £o  I  -—  +  1     +  jl(2  +  5eJ)52  -  2e20]  cos 


2  i  -    2 

+  ^{-9s2  +  8)  cos2y  +  ^l11^  -  6)  cos(</  +  20)  +  -t~(352  -  2)cos(2y  +  20)         (33; 
4  3  24 


+ 


o  c 


2s  K 2 

cos(2y  -  26) -  +  e0 


('2—--  +4-  5s2)  cosy 


+  eo 


<fyl  . 


siny 


c  \    av  )  d62 

In  the  above  equation,  the  cosy  and  siny  terms  would  produce  secular  terms  #siny  and 

6 cosy  in  u\.   The  choice  dy\/d6  =  5s2 /2  —  2  will  eliminate  these  possibilities.   Integrating 

yields 

(34) 


Ifi  =  (  ~Y  ~  2  J  (*  "  tfo)  +  A'3[sin(2y  -  20)  +  sin  2u;0] 

The  term  multiplied  by  A'3  was  added  to  eliminate  secular  terms  in  u2-  The  constant  terms 
in  (34)  were  added  to  satisfy  the  initial  condition  y(0o)  =  60  —  u;0. 
A  solution  to  lowest  order  of  equation  (33)  is  then 


ui«l-T 


+  4  U^-  +  1  j  +  ^[--s2(2  +  5e0)  +  2e2]  cos  2d 


eo 


+  -^(9s2  -  8)  cos2y  +  ^(-l^2  +  6)cos(y  +  20)  +  ;£(-3s2  +  2)  cos(2y  +  26)         (35; 
12  24  24 

+ 


8  c 


2sA'2 

cos(2y  -  2d) +  A'4  cos(y  -  26) 


+  A'5  cos(y  -  0O  +  *>'o)  +  A'6  sin(y  -  60  +  u;0) 

The  term  multiplied  by  K4  was  added  to  eliminate  secular  terms  in  u2-  The  terms  multiplied 
by  A'5  and  J\'e  were  added  to  satisfy  the  initial  conditions  (14)— (16). 

With  all  terms  in  place  to  deal  with  secular  terms,  the  calculations  are  continued  by 
substituting  equations  (28)-(30).  (32).  (34).  and  (35)  into  (26)  and  equating  terms  multiplied 

by  J2: 

sce2(\5s2  -  14)' 


di-i 
~d~6 


A'i  + 


sin(2y-20)  +  . 


(36) 


24(5s2-4) 
We  have  for  brevity  only  indicated  on  the  right  side  of  equation  (36)  the  term  that  would 

produce  secular  terms  in  i2-  Removal  of  this  term  by  making  its  coefficient  zero  determines 

A'].  Equation  (36)  is  then  integrated  to  determine  i2. 

Continuing  the  procedure  by  equating  the  terms  multiplied  by  J2  in  the  expansion  of 

equation  (27)  determines  y2.  A'3.  and  A'4.    Final  values  of  all  the  constants  are  listed  in 

Appendix  I. 
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Knowing  the  solution  for  z(0),  we  can  determine  0,(6)  by  integrating  equation  (8)  and 
applying  the  initial  condition  (18).  The  angle  0.  which  increases  continuously  from  an  initial 
value  #o,  rnay  be  related  to  the  time  t  by  numerically  integrating  (22). 

5     Solution 


Here  we  assemble  the  complete  solution: 

3s2 


5s' 


r    =    p0/|l  +  eocosy-f  J[l-  — +  ejll-— j  +  — (-(2  +  5ejy  +  2eJ)cos2^ 

2  2 

+  ^(9s2  -  8)  cos  2y  +  ^(-1  Is2  +  6)  cos(y  +  20)  +  ^(-3s2  +  2)  cos(2y  +  20) 
12  24  24 

+  -^(3s2-2)cos(2y-20) 
o 

e0[15(2  +  e2)54  -  14(4  +  c20)s2  +  24]  sin  [f  (5s2  -  4)1  sin[0  +  w0] 


12(5s2-4) 
efc2(15s2-  14)  sin  [f  (5s2  -  4)1  sin  [2^0  -  y  (5s2  -  4; 


6(5s2-4) 

2  2    2 

+  ^(3s2  -  2)  cos(y  -  30o  +  3-*)  -  ?SL  cos(j,  -  50o  +  3^-0) 
24  Id 


e252 

7— cos((/  -  0O  +  3^o! 
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(37) 


+  ^(3s2  -  2)  cos(y  -  20o  +  2u;0)  -  ^ 
4  o 


cos(y  -  40o  +  2u>0] 


^(s2  +  l)cos(y  +  2^o)  +  i[(-2  +  5e2)s2  -  2e2]  cos(y  +  0O  +  w0) 
4  o 


1 


+  -[(6  +  5e2)s2  -  4(1  4-  e2)]  cos(y  -  0O  +  **) 

+  ^[-(14  +  5e2)s2  +  2e2]  cos(y  -  30o  +  uq) 

2  2 

+  ^(9s2  -  4)  cos(y  +  30o  -  a*)  +  ^(-7s2  +  6)  cos(y  4-  0O  -  u*>) 
4o  o 


+  77(-5s2  +  4)cos(y-0o-^o) 
lo 


+  — (2s2  -  1 )  cos(y  +  20o)  +  t("3s2  +  ] )  cos^  ~  2^o)  +  t("3s2  +  2) cos  y 


£0 
4 


£0 
4 


+e0s2  cos(0o  +  uq)  +  ~-  cos(30o  -  w0)  +  s2  cos  20o]}  +  p0O{  J2,  J30) 
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y  =  e-uo  +  j(^--2\(6-e0) 


+ 


Je20         f  (-75s6  +  260s4  -  296s2  +  112)  sin  [f  (5s2  -  4)]  cos  [2w0  -  f  (5s2  -  4)] 


24(5s2-4) 


I 


(5s2 -4) 


+J0s2(-15s2  +  14)(15s2-  13)  cos  2^o  I  +  J20J^(15s2-  13)  cos(0o  +  w0) 

2  2 

+  ^-(15s2-13)cos(30o-^o)  +  7r(15s2-  13)cos20o 
6  2 

+  ^[5(9e2  +  34)s4  +  4(9e2  -  34)s2  -  56c2]}  +  0(  J\  J30) 


(38) 


?o  +  scj\  -  cos  20  +  ^  cos(y  +  20) 
2  6 


en                           e2(-15s2  +  14)  sin    ^  (5s2 -4)   sin 
+  -^  cos(y  -  20)  + 


J6  I?.** 


2uJo  _  21  (5s2  -  4) 


12(5s2  -4 


-i  cos  20o  -  j  cos(30o  -  -o)  -  y  cos(0o  +  u>0)  1  +  0(J2,  J30; 


'39: 


n  =  n0  +  cj 


0O  -  0  +  \  sin  20  -  co  sin  y  +  ^  sin(y  +  20)  -  ^  sin(y  -  20)  -  -  sin  20o 

1  0  I  I 


+e0  sin(0o  -  u>0)  -  —  sin(30o  -  u;0)  -  —  sin(0o  +  u?0] 

n  / 


+ 


CJC2        |-2(15s4  -  45s2  +  28)  sin  [f  (5s2  -  4)]  cos  [2^0  -  f  (5s2  -  4) 


12(5s2-4) 


(5s2 -4) 


+70s2(15s2  -  14)  cos  2^o|  +  cJ26  -e0s2  cos(0o  +  w0)  -  ^~  cos(30o  -  ^0; 


-s2  cos  20o  +  ^(7s2  -  4)  +  i(-s2  +  6) 


k  +  o(j\re) 


(40) 


1    /-e     f         r(-3^  +  '■>) 

*    =    U+-       r2    1+J^ —^cos20  +  eo(s2-l) 

ft0  ^0o       l  L  2 
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e0(-4s2  +  3)        ,  „       e0(-2s2  +  l) 

cos  y  +  -^ 1  cos{y  +  20)  +    °l      0         j  cos(y  -  20) 


e^J(15^-14)sin    f  (5s2 -4)   sin   2u,'0  -  <?  (5s2  -  4 


12(5s2-4) 
+s2  -  1  +  -  cos  20o  +  -~-  cos(30o  -  w0)  +  -V  cos(^o  +  u;0) 


(41 


\do  +  ^-0{j2,j2e) 

J  "0 


*0 

In  obtaining  the  equations  (37)— (41 ),  use  has  been  made  of  trigonometric  formulas 
to  simplify  terms  containing  the  factor  5s2  —  4  in  the  denominator.  In  the  form  given, 
these  terms  can  clearly  be  seen  to  approach  a  finite  limit  at  the  ''critical  inclination'* 
i0  =  sin-1  \/4/5  =  63°26'  or  116°34'.  Hence  the  solution  is  actually  valid  for  all  values 
of  io.  If  |?o  —  sin-1  v/4/5|  <  J,  the  formulas  (37)— (41 )  can  still  be  used  by  letting  bs2  —  4  =  J, 
or  the  limiting  forms  for  io  — >  sin-1  v/4/5  can  be  used. 

To  check  the  solution,  we  can  see  if  the  specific  mechanical  energy  (18)  of  the  satellite 
remains  constant.  Substitution  of  the  solution  (36)-(37)  into  equation  (10)  plus  (11)  yields 


GM(l-e2)       GMJ2R2(\  -3  sin2  30)       GM 


T  +  V  = K- "- — — = ^  + O(J') 

2p0  2r£  p0 

The  right  side  is  easily  recognized  as  the  value  of  the  specific  mechanical  energy  at  the  initial 
time  t0. 

As  a  further  check  on  the  solution,  we  can  see  if  it  reduces  to  our  previous  results  for 
equatorial  and  polar  orbits,  obtained  by  completely  separate  derivations  (Danielson  and 
Snider,  1989).  Setting  2o  =  0  and  using  the  independent  variable  o.  measured  from  the  line 
07,  we  find  that  equations  (37)— (41 )  reduce  to  equations  (18)-(22)  of  our  previous  paper. 
Setting  i0  =  tt/2  and  using  the  expansion  cos(y  +  Jk)  %  cos  y  —  Jks'my,  we  find  that 
equations  (37)— (41 )  reduce  to  equations  (38)— (41 )  of  our  previous  paper. 

Comparing  the  terms  in  the  O-symbols,  we  see  that  the  relative  error  in  equation  (41) 
may  be  greater  than  that  of  equations  (37)-(40).  Since  the  underlined  terms  in  equations 
(37)-(40)  are  of  this  same  order  of  magnitude,  we  can  drop  the  underlined  terms  except 
when  (37)— (38)  are  used  to  calculate  r  in  equation  (41).   The  relative  error  of  our  solution 
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will  then  still  be  of  order  (0  -  0o)J2. 

If  we  retain  only  the  two-body  solution,  the  relative  error  terms  will  be  of  the  order 
{6  —  60)J .  Here  the  error  in  our  solution,  as  compared  to  the  exact  solution  of  the  equations 
of  motion,  should  be  of  the  order  J  times  the  error  in  the  two-body  solution  (for  an  Earth 
satellite  J  <  .0015). 

6  Comparison  of  Perturbation,  Two-Body,  Numerical,  and  Mea- 
sured Solutions 

In  this  section  we  compare  the  preceding  perturbation  solution,  the  two-body  solution,  a 
completely  numerical  solution  of  the  differential  equations,  and  actual  measured  satellite 
data;  for  more  comparisons  see  the  thesis  of  Krambeck.  The  difference  between  the  position 
vector  r  determined  by  the  numerical  integration  code  or  measured  data  and  the  position 
vector  rref  calculated  from  our  perturbation  solution  or  the  two-body  solution  is  the  error 
Ar: 

Ar  =  r  -  rref 

If  the  errors  (Ar.  A#,  A?,  Aft)  in  the  orbital  parameters  (r.#.  ?,  Q)  are  small,  we  can  estimate 
Ar  from  equation  (4)  and  the  linear  approximation 

A,  **Ar+£A« +£*  +  £*«  (42, 

Or  39  oi  OYi 

It  is  customary  to  decompose  Ar  into  components  (^i , <52.  ^3)  along  the  moving  triad  (Bi.B2,B3): 

Ar  =  <5]Bi  4-  62B2  +  <53B3 

The  component  8\  is  called  the  radial  error,  62  is  the  down  track  error,  and  <53  is  the  cross 
track  error.  Applying  (42)  to  equation  (4),  and  expressing  the  base  vectors  (bi.b2.b3)  in 
terms  of  (B!,B2.B3),  we  obtain  the  following  approximations: 

6i»Ar.         62  %  r(A0  +  cos?'Afi).         63  «r(sin0At  -  cosflsiniAft)  (43) 
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We  obtained  the  numerical  integration  code  UTOPIA  from  the  Colorado  Center  for 
Astrodynamics  Research  located  on  the  campus  of  the  University  of  Colorado.  The  code 
was  specialized  to  the  differential  equations  used  in  this  paper.  We  compared  the  solutions 
for  an  earth  satellite  with  the  following  initial  conditions: 

r0  =  7,386.18  km 

e0  =  .003991 

0o  =  104.05° 

u;0  =  224.38° 

i0  =  90.03° 

Q0  =  322.63° 

t0  =  0 

These  initial  conditions  represent  an  essentially  polar  orbit  at  an  altitude  of  approximately 
1000  kilometers  and  period  about  l|  hours.  For  this  satellite  the  perturbation  and  numerical 
orbits  match  extremely  well  while  the  two-body  orbit  is  grossly  erroneous.  The  magnitude  of 
the  error  in  r  is  shown  in  Figure  2.  Note  that  the  relative  error  in  our  perturbation  solution 
is  2.8J2(0  —  #0)i  and  that  this  error  is  1.1J  times  the  error  in  the  two-body  solution. 

We  obtained  measured  satellite  data  from  the  First  Satellite  Control  Squadron  located 
at  Falcon  Air  Force  Base.  Colorado.  A  near  earth  satellite  processed  the  following  initial 
conditions: 

r0  =  7, 776.58  km 

e0  =  .0003071 

0O  =  149.14° 

u-'o  =  9.57° 

i0  =  98.81° 
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n0    =    37.10° 
to    =    0000Z  26  July  1990 

Again,  the  perturbation  orbit  is  far  superior  to  the  two-body  orbit.  The  radial,  down  track. 
and  cross  track  errors  (61,62,63)  are  shown  in  Figure  3.  Note  that  although  the  perturbation 
solution  produces  only  a  small  improvement  in  the  radial  error,  this  error  is  negligible  in 
comparison  to  the  down  track  error. 

7     Conclusions 

Our  solution  embodies  the  principles  outlined  in  the  introduction.  The  relative  error  of  our 
solution  is  of  order  (6  —  60)J2,  which  is  a  factor  of  J  times  the  relative  error  of  the  two-body 
solution;  our  solution  loses  its  validity  after  an  angular  change  (6  —  60)  of  order  1/J2,  which 
is  a  factor  of  j  longer  than  the  interval  of  validity  of  the  two-body  solution.  Secondly,  our 
solution  is  in  terms  of  classical  orbital  elements;  no  transformation  to  an  alternative  non- 
physical  set  of  elements  is  required.  Finally,  our  solution  is  free  of  singularities  for  all  values 
of  the  initial  orbital  parameters,  including  elliptic,  parabolic,  and  hyperbolic  orbits. 

Our  formulas  should  agree  closely  with  satellite  orbits  whose  dominant  perturbation  is 
the  planet's  oblateness.  Of  course,  the  effects  of  higher-order  terms  in  these  expansions, 
higher-order  terms  in  the  planet's  potential,  and  of  other  perturbation  forces  may  also  be 
important.  The  formulas  will  have  to  be  amended  to  include  these  additional  effects. 

APPENDIX  I:     Values  of  the  Constants  A'l-A'e 


cse2(-15s2  +  14; 
Aj  — 


24(5s2-4) 


1-  sc        on       sceo       ton  \      sce°       f0  \  .   cseg(15-s2  -  14) 

A  2  =  -—cos  2^o —  cos(30o  -u>0) —  cos(0o  +  ^0)  +      .,..  , --— cos2^0 

2d  I  24(5s''  —  4) 

e§(-75s6  +  260.54  -  296s2  +  112) 
A3  = 


48(5s2-4)2 
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A4  =  e0 


[15(e2  +  2).s4-14(e2  +  4).s2  +  24] 
24(5s2-4) 


K5    =     ^(-952  +  8)cos(2^-2^o)  +  ^t(352-2)cos(4^-2^o) 
12  24 

e  e2 

-(e052  +  A'4)  cos(0o  +  ^o)  +  -^(s2  -  2)  cos(30o  -  u;0)  +  -^(-3s2  +  2)  cos  2u>0 

h  8 

-^[5(2  -  e2)62  +  2e2]  cos  260  +  ^(15e2  +  18)52  -  (eg  +  1) 


2  2 

K6    =    -f^(652-5)sin(2^0-2^o)  +  ^(-3^2  +  l)sin(4^-2u;o) 

D  12 


1 


eo 


+  -[e0(--s2  +  1 )  +  2A4]  sin(0o  +  w0)  +  ^(3s2  -  2)  sin(0o  -  ^o) 


e  e2  1 

+  -^(-7s2  +  2)  sin(30o  -  u;0)  +  -f{s2  +  1 )  sin2u;0  +  -[-(5e2  +  2)s2  +  2eg]  sin  20o 
8  4  0 


APPENDIX  II:      Rigorous  Bounds  on  the  Orbit 


It  follows  from  (10)-(12)  that 


„,      „        1   (drY       r   d2r 

r+V=2A      +5^ 


— (1  —  3  sin   j, 


2r 


4r. 


This  can  be  rewritten  in  the  form 

2" 


rfr 


=  4(7+  \>  +  26U/  + 


GMJ2ff 


;3sin2^-  1) 


from  whence  it  follows  that 


d_ 
dr 


dr 
di 


Integrating  from  r(t0)  to  r(t)  yields 


r>  f  »V  <  2(r  +  vy  +  2GMr  -  >™ML  _  ftg  +  iMMI 


It  follows  that 


0  <  2(T  +  V>2  +  2GMr  -  h\[\  -  ^-] 


(44 
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When  T  +  V  <  0,  the  quadratic  polynomial  on  the  right  side  of  (44)  has  the  roots  (exact 
values  can  be  found  from  the  quadratic  formula) 

Tmin  =   T^-[l  +  0(J2)}  ,  rmax  =   -^-[1  +  0(J2)} 

1  +  e0  1  —  e0 

Hence  a  satellite  having  negative  mechanical  energy  remains  for  all  time  within  the  spherical 
annulus  r^n  <  r  <  r^x.  Since  the  position  vector  is  bounded,  we  can  invoke  the  recurrence 
theorem;  i.e.,  the  satellite  will  come  as  close  as  desired  to  its  initial  position  in  a  sufficiently 
long  period  of  time  (as  shown  by  Poincare).  Furthermore,  we  are  guaranteed  of  the  validity 
of  supressing  secular  terms  to  describe  the  orbit  via  perturbation  analysis. 
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Figure  1:     Orbital  geometry. 


Figure  2:     Comparison  of  perturbation,  two-body,  and  numerical  orbits. 


Figure  3:     Comparison  of  perturbation,  two-body,  and  measured  orbits. 
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